/* taken from i.pca */

#include <stdlib.h>
#include <grass/gmath.h>
#include <grass/gis.h>


static int egcmp (const void *pa, const void *pb);


/*!
 * \fn int eigen (double **M, double **Vectors, double *lambda, int n)
 *
 * \brief Computes eigenvalues (and eigen vectors if desired) for
 * symmetric matices.
 *
 * Computes eigenvalues (and eigen vectors if desired) for symmetric matices.
 *
 * \param M Input matrix
 * \param Vectors eigen output vector matrix
 * \param lambda Output eigenvalues
 * \param n Input matrix dimension
 * \return int
 */

int 
eigen ( 
    double **M,	       /* Input matrix */
    double **Vectors,  /* eigen vector matrix -output */
    double *lambda,    /* Output eigenvalues */
    int n              /* Input matrix dimension */
)
{
	int   i,j;
	double **a,*e;

	a = G_alloc_matrix(n,n);
	e = G_alloc_vector(n);

	for(i=0; i<n; i++)
	for(j=0; j<n; j++) 
	  a[i][j] = M[i][j];

	G_tred2(a,n,lambda,e);
	G_tqli(lambda,e,n,a);
	
	/* Returns eigenvectors	*/
	for(i=0; i<n; i++) 
	for(j=0; j<n; j++) 
	  Vectors[i][j] = a[i][j]; 

	G_free_matrix(a);
	G_free_vector(e);

    return 0;
}


/*!
 * \fn int egvorder2 (double *d, double **z, long bands)
 *
 * \brief
 *
 * Returns 0.
 *
 * \param d
 * \param z
 * \param bands
 * \return int
 */

int 
egvorder2(double *d, double **z, long bands)
{
  double *buff;
  double **tmp;
  int i, j;

  /* allocate temporary matrix */
  buff = (double *) G_malloc(bands * (bands + 1) * sizeof(double));
  tmp = (double **) G_malloc(bands * sizeof(double *));
  for (i = 0; i < bands; i++)
    tmp[i] = &buff[i * (bands + 1)];

  /* concatenate (vertically) z and d into tmp */
  for (i = 0; i < bands; i++) {
    for (j = 0; j < bands; j++)
      tmp[i][j+1] = z[j][i];
    tmp[i][0] = d[i];
  }

  /* sort the combined matrix */
  qsort(tmp, bands, sizeof(double *), egcmp);

  /* split tmp into z and d */
  for (i = 0; i < bands; i++) {
    for (j = 0; j < bands; j++)
      z[j][i] = tmp[i][j+1];
    d[i] = tmp[i][0];
  }

  /* free temporary matrix */
  G_free(tmp);
  G_free(buff);

  return 0;
}

  
/*!
 * \fn int transpose2 (double **eigmat, long bands)
 *
 * \brief
 *
 * Returns 0.
 *
 * \param eigmat
 * \param bands
 * \return int
 */

int 
transpose2(double **eigmat, long bands)
{
  int i, j;

  for (i = 0 ; i < bands ; i++)
    for (j = 0 ; j < i; j++)
    {
      double tmp = eigmat[i][j];

      eigmat[i][j] = eigmat[j][i];
      eigmat[j][i] = tmp;
    }

    return 0;
}


static int
egcmp (const void *pa, const void *pb)
{
  const double *a = *(const double * const *)pa;
  const double *b = *(const double * const *)pb;

  if (*a > *b) return -1;
  if (*a < *b) return 1;

  return 0;
}
